This script assembles various environmental layers into a common 30m grid for the Cape Peninsula focussing only on the Working for Water plots. It also calculates veg age based on the fire data.
Import raster of an index grid (ig) to spatially connect all the datasets.
ig=raster(paste0(datadir,"clean/indexgrid_landsat_30m.grd"))
rv=readOGR(dsn=paste0(datadir,"raw/VegLayers/Vegetation_Indigenous_Remnants"), layer="Vegetation_Indigenous_Remnants")
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/glennmoncrieff/Documents/GIS/Postfire_test/data/raw/VegLayers/Vegetation_Indigenous_Remnants", layer: "Vegetation_Indigenous_Remnants"
## with 3428 features
## It has 7 fields
#remnant veg layer - readOGR() reads shapefiles
#rv; summary(rv$National_); summary(rv$Subtype); summary(rv$Community); levels(rv@data$National_)
rv_meta=data.frame(1:length(levels(rv@data$National_)), levels(rv@data$National_)) #save VegType metadata
colnames(rv_meta)=c("ID", "code") #rename columns
write.csv(rv_meta, "data/vegtypecodes.csv", row.names=F)
# reproject to the CRS of the Landsat index grid (UTM 34S)
rv=spTransform(rv,CRS(proj4string(ig)))
Extract the national veg types from the veg layer into a 30m raster based on the index grid
rvrfile="data/vegtypes_landsat_30m.tif"
## note the if(!file.exists)) first checks if the file already exists so you don't rerun this everytime you run the script...
if(!file.exists(rvrfile))
rvr=rasterize(rv, ig, field=c("National_"), fun="max",file=rvrfile) #get national veg type for each cell
## read it back in and 'factorize' it
rvr=raster(rvrfile)
rvr=as.factor(rvr)
rv_meta$code=as.character(rv_meta$code)
levels(rvr)=rv_meta[levels(rvr)[[1]]$ID,]
levelplot(rvr,col.regions=rainbow(nrow(rv_meta),start=.3))
Count number of veg types for each cell (i.e. ID mixed cells)
rvcfile="data/count_vegtypes_landsat_30m.tif"
if(!file.exists(rvcfile))
rvc=rasterize(rv, ig, field=c("National_"), fun="count",file=rvcfile)
rvc=raster(rvcfile)
Are there any mixed cells?
table(values(rvc))
##
## 1 2
## 337337 5617
alien_plots<-read.csv(paste0(datadir,"clean/alienplots.csv"),stringsAsFactors = FALSE)
alien_xy<-cbind(as.numeric(alien_plots$x),as.numeric(alien_plots$y))
alien_xy<-SpatialPoints(alien_xy,crs(ig))
alien_data<-SpatialPointsDataFrame(alien_xy,alien_plots[,4:8])
##load plot clearing data
alien_n <- readOGR(dsn=paste0(datadir,"raw/AlienPlots/Nbal Analysis_TMNP_20160314"),"Nbal Analysis_TMNP_North_20160314")
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/glennmoncrieff/Documents/GIS/Postfire_test/data/raw/AlienPlots/Nbal Analysis_TMNP_20160314", layer: "Nbal Analysis_TMNP_North_20160314"
## with 306 features
## It has 123 fields
alien_c <- readOGR(dsn=paste0(datadir,"raw/AlienPlots/Nbal Analysis_TMNP_20160314"),"Nbal Analysis_TMNP_Central_20160314")
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/glennmoncrieff/Documents/GIS/Postfire_test/data/raw/AlienPlots/Nbal Analysis_TMNP_20160314", layer: "Nbal Analysis_TMNP_Central_20160314"
## with 436 features
## It has 123 fields
alien_s <- readOGR(dsn=paste0(datadir,"raw/AlienPlots/Nbal Analysis_TMNP_20160314"),"Nbal Analysis_TMNP_South_20160314")
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/glennmoncrieff/Documents/GIS/Postfire_test/data/raw/AlienPlots/Nbal Analysis_TMNP_20160314", layer: "Nbal Analysis_TMNP_South_20160314"
## with 206 features
## It has 123 fields
alien_n<-spTransform(alien_n,CRS(proj4string(ig)))
alien_c<-spTransform(alien_c,CRS(proj4string(ig)))
alien_s<-spTransform(alien_s,CRS(proj4string(ig)))
nn<-seq(1,nrow(alien_n@data),by=1)
nc<-seq(max(nn)+1,max(nn)+nrow(alien_c@data),by=1)
ns<-seq(max(nc)+1,max(nc)+nrow(alien_s@data),by=1)
nn1 <- spChFIDs(alien_n, as.character(nn))
nc1 <- spChFIDs(alien_c, as.character(nc))
ns1 <- spChFIDs(alien_s, as.character(ns))
temp<-spRbind(nn1,nc1)
vegclear<-spRbind(temp,ns1)
fi=readOGR(dsn=paste0(datadir,"raw/Fire"), layer="CapePenFires") #Cape Peninsula fires history layers 1962-2007
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/glennmoncrieff/Documents/GIS/Postfire_test/data/raw/Fire", layer: "CapePenFires"
## with 4578 features
## It has 9 fields
fi=spTransform(fi,CRS(proj4string(ig)))
### Extract fire history data and convert to a 30m raster
fi$STARTDATE[which(fi$STARTDATE==196201001)]=19620101#fix an anomalous date...
#Raster showing total numbers of fires in each grid cell
ficfile="data/fires_number_1962to2007_landsat_30m.tif"
if(!file.exists(ficfile))
fic=rasterize(fi, ig, field=c("YEAR"), fun="count",file=ficfile)
fic=raster(ficfile)
## Set up a "Days since 1 January 1960" column on the fire polygons
sdate=fi$STARTDATE #get the unique list of years in which fires occurred
sdate[which(substr(sdate, 5, 8)=="0000")]=sdate[which(substr(sdate, 5, 8)=="0000")]+1231 #set those fires to occur on the 31st Dec - need to check these and fix
ddate=as.Date(as.character(sdate), format="%Y%m%d")
sdate=as.numeric(ddate-as.Date("19600101", format="%Y%m%d"))
fi$Day=sdate
fi$DayDate=ddate
sdate=sort(unique(sdate))
getNDVI=function(file,datefile,prefix){
ndvi=stack(paste0(datadir,"raw/NDVI/",file))
NAvalue(ndvi)=0
offs(ndvi)=-2
gain(ndvi)=.001
dates=as.Date(read.csv(paste0(datadir,"raw/NDVI/",datefile),header=F)[,1])
names(ndvi)=paste0(prefix,sub("-","",dates))
ndvi=setZ(ndvi,dates)
}
Now use the function to read in the data and add the relevant metadata.
l4=getNDVI(file="2016045_v1g_LT4_L1T_TOA_daily__1988-03-08-1992-11-14.tif",
datefile="LT4_L1T_TOA2.csv",prefix="L4_")
l5=getNDVI(file="2016045_v1g_LT5_L1T_TOA_daily__1984-06-09-2011-04-17.tif",
datefile="LT5_L1T_TOA2.csv",prefix="L5_")
l7=getNDVI(file="2016045_v1g_LE7_L1T_TOA_daily__1999-08-30-2016-03-21.tif",
datefile="LE7_L1T_TOA2.csv",prefix="L7_")
l8=getNDVI(file="2016045_v1g_LC8_L1T_TOA_daily__2013-05-24-2016-03-13.tif",
datefile="LC8_L1T_TOA2.csv",prefix="L8_")
Let’s check out one of the LANDSAT objects. Raster provides a summary by just typing the object’s name:
l7
## class : RasterStack
## dimensions : 1929, 674, 1300146, 229 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 249990, 270210, 6189390, 6247260 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=34 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : L7_199908.30, L7_199910.01, L7_199911.02, L7_199912.04, L7_200001.21, L7_200004.10, L7_200006.13, L7_200007.31, L7_200010.03, L7_200010.19, L7_200011.20, L7_200102.24, L7_200103.12, L7_200103.28, L7_200104.29, ...
## min values : -34.768, -34.768, -34.768, -34.768, -34.768, -34.768, -34.768, -34.768, -34.768, -34.768, -34.768, -34.768, -34.768, -34.768, -34.768, ...
## max values : 30.767, 30.767, 30.767, 30.767, 30.767, 30.767, 30.767, 30.767, 30.767, 30.767, 30.767, 30.767, 30.767, 30.767, 30.767, ...
## time : 1999-08-30 - 2016-03-21 (range)
And a plot of a few different dates:
yearind=which(getZ(l7)%in%getZ(l7)[1:5])
levelplot(l7[[yearind]],col.regions=cndvi()$col,cuts=length(cndvi()$at),at=cndvi()$at,layout=c(length(yearind),1),scales=list(draw=F),maxpixels=1e5)
There is some temporal overlap between sensors, let’s look at that:
tl=melt(list(l4=getZ(l4),l5=getZ(l5),l7=getZ(l7),l8=getZ(l8)))
xyplot(as.factor(L1)~value,data=tl,pch=16,groups=as.factor(L1),asp=.15,lwd=5,ylab="LANDSAT Satellite",xlab="Date")
Timeline of LANDSAT data by sensor
There are several ways these data could be combined.
The individual scenes could be assessed for quality (cloud contamination, etc.), sensors could be weighted by sensor quality (newer=better?).
Today, we’ll simply combine (stack) all the available observations for each pixel.
ndates=unique(tl$value)
## write that to disk for later use
write.table(ndates,file="data/ndates.csv",row.names=F,col.names=F,sep=",")
### concatenate all sensors to a single raster stack (this just links the filenames)
undvi=stack(l4,l5,l7,l8)
undvi=setZ(undvi,c(getZ(l4),getZ(l5),getZ(l7),getZ(l8)))
### Sort them by date (z)
ndvi=undvi[[order(getZ(undvi))]]
ndvifile="data/ndvi_wfw.Rdata"
if(!file.exists(ndvifile)){
save(ndvi,file=ndvifile)
# ## create a new ndvi layer
# ndvifile="data/ndvi_landsat_30m.tif"
# if(!file.exists(ndvifile)){
# writeRaster(ndvi,filename=ndvifile)
# }
# ## load it
# ndvi=raster(ndvifile)
}
load(ndvifile)
yearind=which(getZ(ndvi)%in%getZ(l7)[1:5])
levelplot(ndvi[[yearind]],col.regions=cndvi()$col,cuts=length(cndvi()$at),
at=cndvi()$at,margin=F,scales=list(draw=F),
names.attr=getZ(ndvi)[yearind],maxpixels=1e4)
Merged annual maximum LANDSAT NDVI
Here we will define the subset of cells that we will explore further. In this case it is the WfW plots. You can fiddle with these settings to include fewer (or more) cells. If your computer is slow, you may want to subset this further.
## load data for masking
cover=raster(paste0(datadir,"clean/landcover2009_landsat_30m.gri"))
maskfile="data/mask_landsat_30m.tif"
if(!file.exists(maskfile)){
mask=overlay(cover,fic,fun=function(x,y) x==1&y>0,filename=maskfile)
}
mask=raster(maskfile)
## load additional covariate data
tmax=raster(paste0(datadir,"clean/Tmax_jan_mean.gri"))
tmin=raster(paste0(datadir,"clean/Tmin_jul_mean.gri"))
tpi=raster(paste0(datadir,"clean/tpi500.gri"))
dem=raster(paste0(datadir,"clean/dem_landsat_30m.gri"))
janrad=raster(paste0(datadir,"clean/janrad.gri"))
julrad=raster(paste0(datadir,"clean/julrad.gri"))
aspect=raster(paste0(datadir,"clean/aspect.gri"))
### Make a dataframe of all spatial data
## Beware, this approach will only work if your data are all in identical projection/grid/etc.
maskids=extract(mask,alien_xy,cellnumbers=TRUE)
maskids=maskids[,1]
cleardat=data.frame(
id=extract(ig, maskids),
coordinates(ig)[maskids,],
extract(vegclear,alien_xy) #maskids wont work here becasue vegclear is not a raster
)
sdat=data.frame(
id=extract(ig, maskids),
coordinates(ig)[maskids,],
veg=extract(rvr, maskids),
cover=extract(cover, maskids),
tmax=extract(tmax, maskids),
tmin=extract(tmin, maskids),
janrad=extract(janrad, maskids),
julrad=extract(julrad, maskids),
aspect=extract(aspect, maskids),
dem=extract(dem, maskids),
tpi=extract(tpi, maskids)
# firecount=extract(fic, maskids)
)
sdat<-cbind(sdat,alien_data@data)
kable(head(sdat))
| id | x | y | veg | cover | tmax | tmin | janrad | julrad | aspect | dem | tpi | Site | Altitude | Geology | Vegetation | Aliens |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 510630 | 262335 | 6224535 | 16 | 1 | 26.75148 | 8.804244 | 8564.083 | 3034.583 | 4.348300 | 295.9442 | -5.746557 | 1 | Med | TMS | Fynbos | No Aliens |
| 503866 | 261615 | 6224835 | 16 | 1 | 27.39490 | 8.695924 | 8577.000 | 2056.750 | 2.866099 | 310.1948 | 3.405393 | 2 | Med | TMS | Fynbos | No Aliens |
| 509959 | 262425 | 6224565 | 16 | 1 | 27.25568 | 9.311305 | 8444.250 | 3097.167 | 4.679824 | 326.9142 | 34.642132 | 3 | Med | TMS | Fynbos | No Aliens |
| 488295 | 259545 | 6225525 | 16 | 1 | 24.47923 | 7.349550 | 8333.167 | 2424.000 | 2.976255 | 535.4271 | 64.628049 | 4 | High | TMS | Fynbos | Invaded |
| 508605 | 262245 | 6224625 | 16 | 1 | 27.00466 | 8.638811 | 8858.917 | 3093.083 | 4.129409 | 286.9676 | -14.465683 | 5 | Med | TMS | Fynbos | No Aliens |
| 481553 | 259485 | 6225825 | 16 | 1 | 25.35269 | 7.169548 | 8704.250 | 2204.333 | 2.700882 | 560.4637 | 39.230657 | 6 | High | TMS | Fynbos | Cleared |
It’s often easier to work with data in ‘long’ format where there is one row for each observation and another column indicating what the observation is. Let’s melt the data to ‘long’ format.
#extract ndvi data
ftdatw="data/tdatw.Rdata"
if(!file.exists(ftdatw)){
tdatw=data.frame(
id=extract(ig, maskids),
# extract(age, maskids),
extract(ndvi,maskids)
)
save(tdatw,file=ftdatw)
}
load(ftdatw)
kable(tdatw[1:10,1:10])
| id | L5_198406.09 | L5_198409.29 | L5_198410.15 | L5_198504.09 | L5_198611.06 | L5_198611.22 | L5_198612.24 | L5_198701.09 | L5_198702.26 |
|---|---|---|---|---|---|---|---|---|---|
| 510630 | 0.457 | 0.116 | -2.000 | 0.377 | -2 | -2 | 0.423 | 0.457 | -2 |
| 503866 | 0.439 | 0.103 | -2.000 | 0.358 | -2 | -2 | 0.341 | 0.368 | -2 |
| 509959 | 0.339 | 0.078 | -2.000 | 0.126 | -2 | -2 | 0.222 | 0.148 | -2 |
| 488295 | 0.682 | 0.074 | -2.000 | 0.529 | -2 | -2 | -2.000 | 0.450 | -2 |
| 508605 | 0.546 | 0.103 | 0.385 | 0.342 | -2 | -2 | 0.397 | 0.395 | -2 |
| 481553 | 0.509 | 0.055 | -2.000 | 0.348 | -2 | -2 | -2.000 | 0.390 | -2 |
| 521371 | 0.649 | 0.203 | 0.535 | 0.596 | -2 | -2 | 0.631 | 0.642 | -2 |
| 532137 | 0.629 | -2.000 | 0.322 | 0.383 | -2 | -2 | 0.364 | 0.444 | -2 |
| 530791 | 0.610 | -2.000 | 0.566 | 0.599 | -2 | -2 | 0.527 | 0.538 | -2 |
| 475386 | 0.523 | -2.000 | -2.000 | 0.523 | -2 | -2 | 0.377 | 0.391 | -2 |
tdatl=melt(tdatw,id.var="id")
tdatln=cbind.data.frame(lab=levels(tdatl$variable),
do.call(rbind,strsplit(as.character(levels(tdatl$variable)),"_")))
tdatln$type=ifelse(tdatln[,"1"]%in%c("L4","L5","L7","L8"),"ndvi","age")
tdatl[,c("type","date")]=tdatln[match(tdatl$variable,tdatln$lab),4:3]
tdatl$miss=paste(substr(as.character(tdatl$variable), 1, 2))
tdat=dcast(tdatl,id+date+miss~type,value.var="value")
## convert date to proper format
n = 5 #where to insert separator
tdat$date=paste(substr(as.character(tdat$date), 1, 5-1), ".", substr(as.character(tdat$date), n, nchar(as.character(tdat$date))), sep = "")
## convert year from a factor to numeric
tdat$date=as.Date(as.character(tdat$date),"%Y.%m.%d")
kable(head(tdat))
| id | date | miss | ndvi |
|---|---|---|---|
| 319827 | 1984-06-09 | L5 | 0.761 |
| 319827 | 1984-09-29 | L5 | 0.269 |
| 319827 | 1984-10-15 | L5 | -2.000 |
| 319827 | 1985-04-09 | L5 | 0.596 |
| 319827 | 1986-11-06 | L5 | 0.498 |
| 319827 | 1986-11-22 | L5 | -2.000 |
#extract fire data for each wfw plot
fi_plot <- extract(fi,alien_xy)
#create a seq of date from earliest NDVI to most recent
alldate <- seq.Date(min(tdat$date),max(tdat$date),by="day")
firefile="data/fire.age_wfw.Rdata"
if(!file.exists(firefile)){
fire.age <- matrix(,length(alien_xy),length(alldate))
#loop through plots
for (i in 1:length(alien_xy)){
# get fire dates for each pixel
dates.b <- filter(fi_plot, point.ID == i) %>% select(DayDate)
#loop through days - ouch
for (j in 1:length(alldate)){
if(any(!is.na(dates.b)) && ((alldate[j]-min(dates.b[,1],na.rm=T))>=0)){
fire.diff <- alldate[j] - dates.b[,1]
fire.diff.pos <- fire.diff[which(fire.diff>=0)]
fire.age[i,j] <- min(fire.diff.pos)
}
}
}
fire.age<-as.data.frame(fire.age)
names(fire.age) <- alldate
save(fire.age,file=firefile)
}
load(firefile)
#create data frame
fi.plot.dat <- data.frame(
id=extract(ig, alien_xy),
fire.age)
names(fi.plot.dat) <- c("id",as.character(alldate))
fi.age<-melt(fi.plot.dat,id.vars="id")
fi.age$variable<-as.Date(fi.age$variable)
#join NDVI and age data
tdat.all<-left_join(tdat, fi.age, by = c("id" = "id", "date" = "variable"))
tdat.all.com<-tdat.all[complete.cases(tdat.all),]
## check it out
kable(head(tdat.all),row.names = F)
| id | date | miss | ndvi | value |
|---|---|---|---|---|
| 319827 | 1984-06-09 | L5 | 0.761 | NA |
| 319827 | 1984-09-29 | L5 | 0.269 | NA |
| 319827 | 1984-10-15 | L5 | -2.000 | NA |
| 319827 | 1985-04-09 | L5 | 0.596 | NA |
| 319827 | 1986-11-06 | L5 | 0.498 | NA |
| 319827 | 1986-11-22 | L5 | -2.000 | NA |
#join clearing data
clearsmall <- select(cleardat, id,point.ID,First_Date, First_End_, Second_End, Third_End_, Fourth_End, Fifth_End_)
clearsmall$First_Date <- as.Date(clearsmall$First_Date, format="%Y%m%d")
clearsmall$First_End_ <- as.Date(clearsmall$First_End_, format="%Y%m%d")
clearsmall$Second_End <- as.Date(clearsmall$Second_End, format="%Y%m%d")
clearsmall$Third_End_ <- as.Date(clearsmall$Third_End_, format="%Y%m%d")
clearsmall$Fourth_End <- as.Date(clearsmall$Fourth_End, format="%Y%m%d")
clearsmall$Fifth_End_ <- as.Date(clearsmall$Fifth_End_, format="%Y%m%d")
#matrix with a row for every plot and columns up to max fires
fmat<-matrix(,nrow=length(unique(fi_plot$point.ID)),ncol= max(count(fi_plot,point.ID)[,2]))
fnum<-numeric(length(unique(fi_plot$point.ID)))
for(i in 1:length(unique(fi_plot$point.ID))){
temp <- filter(fi_plot, point.ID == i) %>% select(DayDate)
temp <- as.character(temp[,1])
fmat[i,1:length(temp)]<-temp
fnum[i] <- i
}
fmat<-as.data.frame(fmat)
fmat <- cbind(fnum,fmat)
names(fmat) <- c("point.ID","fire1","fire2","fire3","fire4")
fmat$fire1 <- as.Date(fmat$fire1, format="%Y-%m-%d")
fmat$fire2 <- as.Date(fmat$fire2, format="%Y-%m-%d")
fmat$fire3 <- as.Date(fmat$fire3, format="%Y-%m-%d")
fmat$fire4 <- as.Date(fmat$fire4, format="%Y-%m-%d")
clear.all.dat <- left_join(clearsmall,fmat,by="point.ID")
read vertical lines show clearing, purple vertical lines show fires
read vertical lines show clearing, purple vertical lines show fires
showing only L5
showing only L5
showing only L5
showing only L5
## drop the 'wide' version
save(clear.all.dat,tdat.all,file="data/modeldata_nofire_wfw.Rdata")